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Abstract 

O ■ We propose a new Monte-Carlo method for calculation of the Casimir forces. Our 

^ ! method is based on the formalism of noncompact lattice quantum electrodynamics. 

^1 This approach has been tested in the simplest case of two ideal conducting planes. 

After this the method has been applied to the calculation of the lateral Casimir forces 
between two ideal conducting rectangular gratings. We compare our calculations 
with the results of PFA and "Optimal" PFA methods. 
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Introduction 

^ I The Casimir effect attracts significant attention in the last few years due to its important role 
V/-^ ■ in micro- and nano- mechanics [IJ. There are two main problems in this area: calculation of 
O ■ Casimir forces for bodies that are interesting from experimental point of view; consideration of 
. both the complicated shape and electromagnetic properties of the interacting bodies. It is also 
^ ' important to take into account temperature effects. Radiative corrections are small in QED 
[2] . But they may play more important role in case of non-abelian fields. Our approach makes 
it possible to take them into account. 
' There are many numerical methods have been proposed for calculation of the Casimir energy. 

■ Let us discuss some of them. 
- - -' One of the most efficient methods for Casimir energy calculation was derived in [3]. In this 

case the stress tensor and the net force on a body is obtained from Euclidean Green's function. 
It can be calculated numerically by using of standard electromagnetic methods (for example, 
by Finite Elements Method or Boundary Elements Method). 

Another interesting numerical method for the Casimir energy calculation was proposed in 
[3]. This worldline algorithm is formulated for calculations of Casimir forces induced by scalar- 
field fluctuations with Dirichlet boundary conditions for various geometries. Unfortunately this 
method was not extended to Neumann boundary conditions. 

Our approach is based on Quantum Field Theory on the lattice. It is one of the common 
frameworks in modern science. It is based on the Monte-Carlo calculations and provides an 
efficient tool to study all the above-mentioned problems. We have already studied the Casimir 
interaction of the Chern-Simons surfaces by using lattice QED O [6]. In the present paper 
we consider Monte-Carlo calculation of the Casimir interaction between ideal conductors and 
between dielectric bodies. As an example, we study the Casimir effect for two rectangular 
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gratings. The normal part of the vacuum force in this system has already been calculated [?]. 
But there exist also lateral forces, which are of significant experimental and technical interest. 

1 Ideal conductors and dielectric bodies in QED on the 
lattice 

1.1 Conductor boundary conditions 

In this paper we use the 4-dimensional hypercubical lattice and the simplest action of noncom- 
pact lattice QED in Euclidean time: 

X fl<U 

where the link and plaquette variables are defined as follows: 

U^{x) = eaA^, (2) 

9^,{x) = A^U,{x) - A,f/^(x), 

Af,Uu{x) = Uu{x + jj) - Uu{x), 

where a is the lattice step, (3 = 1/e^. The lattice is formed by the sites. They are vertices of 
4-dimensional cubes. The links are edges of these cubes and plaquettes are their sides. The 
functional integrals for physical quantities are calculated by means of the Monte-Carlo method. 
It means that we generate field configurations (sets of link variables) with statistical weight 

. After it the functional integrals are calculated as field configurations averages. 

Our first task is formulation of the ideal conductor boundary condition on the lattice. Both 
electric and magnetic fields are pushed away from the interior of conductor. So the boundary 
conditions for the fields can be written as: 

E||U = 0, H^\s = 0. (3) 

Because the hypercubical lattice is used in our approach, we will approximate a surface by a 
set of plane polygons. Each polygon consists of plaquettes. 

Let us consider a surface that consists of one polygon (normal to the z-axis). More com- 
plicated surfaces can be treated similarly. Expressions ([3]) are the boundary conditions for the 
field strength tensor F^^. Due to the connection between lattice variables and fields ([2]), the 
boundary conditions can be reformulated in terms of plaquette variables: 

e^lix) = 0, ^42 (X) = 0, ^12 (X) = . (4) 

All plaquettes in these conditions form 3-dimensional sublattice that describes the conductor 
surface. Pure gauge is the only type of field configuration that satisfies this conditions. It means 
that any link variable inside this sublattice can be written as: 

Ui = a{x + i) — a{x), (5) 

where a{x) is an arbitrary function defined on the sites of 3-dimensional sublattice. This scheme 
is valid for any surface. 

The functional integrals should be rewritten, taking into account the conductor boundary 
conditions. Our functional integrals must be formulated in terms of independent field variables 
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Figure 1: Two parallel conducting planes: sublattices of links with zero value link variables. 

A^{x), but in our case some of them are not independent due to the conditions ([3]). This problem 
can be solved by choosing the function a{x) as an independent variable which parameterizes 
the fields on the boundary surface. 

So the functional integral for the partition function takes the form: 



observable can be rewritten in terms of lattice variables: SEuci.Lat. and F\U^. Then 

the functional integral is Boltzmann-type field configuration average with statistical weight 
exp(— 5). The lattice field configurations consist of the link variables [/^ and the site variables 
a{x). The link variables describe the electro-magnetic field outside the boundary surface. The 
site variables are defined on the lattice sites that belongs to boundary surface. The function 
a{x) parameterizes the pure gauge fields on this surface. 

The pure gauge field on the surface can be transformed into the zero- value fields a{x) = 
by means of a suitable lattice gauge transformation. After such transformation, the action 
SEuci.Lat. and the observable F do not change their values. It means that in order to implement 
the ideal conductor boundary conditions it is sufficient to generate field configurations with 
zero-value link variables on the boundary surface. 

Let us consider the following example: the Casimir interaction of two parallel conducting 
planes. 

The boundary links are always taken equal to zero during the generation procedure. In 
Fig. 1, these boundary links are bold-marked. The boundary links form two 3-dimensional 
sublattices a.t z = and z = N, where is lattice size. The 3-dimensional structure of these 
sublattices is shown in Fig. 1 as well. These boundary conditions can be written explicitly as 



Z 




where A is the vector potential of the fields outside the boundary surface. 

Physical quantities are calculated by means of the Euclidean functional integral: 




Ui{x,y,z,t) =0 if i = 1,2,4 and z = 0,N. 
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The periodic boundary conditions are used for the other three coordinates {x, y and t). 

This approach can be used to implement the ideal conductor boundary conditions for any 
surface shape. 

1.2 Description of dielectric bodies in QED on the lattice 

The simplest dielectric (we neglect anisotropy and dielectric constant frequency dependence) 
can be described by the following action: 

3 

V V 

Here V is the 4-dimensional volume occupied by dielectric body and e is the static dielectric 
constant. This action can be rewritten in terms of lattice quantities: 

^ = f f E E ^^^(^) + E f ^ E ^o^(^) + E ^J(^) 

\x<^V M<'' x&V \ i=l i<j 

We shall use it for field configurations generation. The dielectric constant frequency dependence 
is very important for real materials and has to be taken into account. Firstly, we fix the gauge: 

A4 = . 

Then we make the Fourier transform: 

Ai = ^= / Aie-^^dw, 1,2,3. 
v27r J 

Here u is the imaginary frequency, since we make the Fourier transform in Euclidean time. 
After this, the action becomes "block-diagonal": 

-^eucl = IJ drduj (f^^ + F^, + F^, + e{tu)u\Al + Al + Al)) . 

So we can generate field configurations for each ou independently and take into account the e{uj) 
dependence. 

2 The lattice observable for the ground state energy 

The corresponding lattice observable should give (after the Monte-Carlo averaging) the ground 
state energy of a system: 

/ (0|T°0|0)(if . (7) 
Jv 

Our task now is to define such lattice observable. 

Lattice system can be treated as a certain A^-dimensional quantum-mechanical system where 
the number of degrees of freedom N tends to infinity. For simplicity we will firstly consider the 
one-dimensional problem of quantum particle moving in potential V{x). The lattice formulation 
of this task: the path in functional integral is a set of the coordinate values at different moments 
of time {x{ti) — Xi, i — 0...N). The boundary conditions in the Euclidean time direction are 
taken periodic: xq — x^. Such a path is an analog of the Euclidean field configuration in 
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Lattice Field Theory. The Monte-Carlo method for functional integral calculation is based 
on generation of a set of such lattice configurations with statistical weight proportional to 
exp(~'^™'='t'^'l/^*). The Euclidean action is 

h V{Xr. 




'-'eucl 

where a is the time step. It is well known [10] that any coordinate x„ is distributed with 
the probability density of the vacuum state. It means that if one calculates the configuration 
average of a certain lattice observable (for example, the potential energy), one gets the vacuum 
average of this observable (the average potential energy): 



{0\V\0) = {V{x^)) = 



Nconf 

' con J 

In order to obtain the full energy of the system, we need to calculate the average of the 
kinetic energy as well. This problem was studied in [H], the corresponding lattice observable 
for the kinetic energy is given by 



(oirio) = (-- ^^"+^ , + —). 



The naive expression for the kinetic energy lattice observable (the first part in the sum ([8])) 
diverges in the continuous limit due to fractal structure of the trajectories. This divergence 
is directly cancelled in ([8]). This procedure can be performed for any lattice models (and for 
lattice field theory too), and thus we have the method for direct calculation of the full energy 
on the lattice. 

Let us consider noncompact lattice QED with the Hamiltonian density: 

\i=l i<j / 

where the "field momentum" tTj = Fq^ is the conjugated quantity to the field Ai. The vacuum 
expectation value of the second part of the expression ([9]) can be calculated directly by the 
field configuration averaging. The averaging of the first part of (Q is performed by the same 
way as the kinetic part of one- dimensional quantum mechanical system: we calculate the same 
observable, but with opposite sign at 7rj(x)^. The final expression for vacuum expectation value 
of the full energy reads 

(10) 



The renormalization procedure for the lattice observable f lTOj) has three parts. Firstly, the 
singular contribution in (|8]) that corresponds to fractal structure of the trajectories has to be 
eliminated. So we calculate the energy density (0|T''°(x) |0) on the free lattice and then subtract 
it from the energy density calculated in presence of the boundaries. 



(0|T""(f)|0)bound - (0|T""(x)|0) ^ee lattice- 

As a result, one gets the difference between the Hamiltonian densities with and without inter- 
acting bodies. It is important to note that (0|T''°(f)|0)fj.gg lattice J^^^ ^ constant. 
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But this is not the end of the story. Any bodies have the Casimir self-energies which are 
connected with their boundaries. The Casimir self-energy calculation is an important task, 
e.g., for the dynamical Casimir effect, but in our case this type of the Casimir energy should 
be eliminated because we are interested only in interaction between bodies. 

Finally, we have the following three stages of the renormalization procedure: 

1. Calculation of the energy density (0|T°°(x)|0) in presence of the boundary conditions (in- 
teracting bodies) and subtraction of the "fractal divergence" constant (0|T°°(^)|0) fj.gg lattice- 

2. Calculation of the energy density (0|T°°(x)|0)ggif_gj^gj.gy for all interacting bodies individ- 
ually. 

3. Subtraction of this Casimir self-energy density (0|T°°(a;)|0)ggif_gj^gPgy from the Casimir 
energy density of the interacting bodies after elimination of the "fractal divergence". 
As a result, one obtains the renormalized Casimir energy densuty for the system of the 
interacting bodies. 

3 Results of the numerical calculations 

We use the standard "heat bath" method [TH] for field configurations generation. Also we use 
the anisotropic lattice to improve the accuracy of our calculations. The anisotropic lattice 
formalism is similar to the one applied in [11] to study the finite temperature effects, but in 
the present work we deform 2;-direction instead of t-direction. 

All link variables are connected with the field components in the usual way: 

f/^ = eaA^,/i = 1,2,4 

except links in 2;-direction: 

U3 = aea A3. 

Here a is the lattice step in all directions except z; aa is the lattice step in 2;-direction. The 
action can be written as 

^ = f{"E E ^;.(-) + ^EE^3.(a:)}. (11) 

The continuous limit for the new observable is different from the traditional methods used, 
for example, in lattice QCD. It is well known that in the continuous theory the Casimir energy 
(without radiative corrections) is independent of the electron charge. This property is exactly 
reproduced in noncompact lattice electrodynamics for the following reasons. Firstly, the phase 
transition is absent in this theory, thus /3 in ( ITT]) plays only the role of parameter that deter- 
mines the numerical values of link variables. Therefore, if one calculates the Casimir energy in 
noncompact electrodynamics without fermions, the value of /3 can be chosen arbitrarily within 
sufficiently broad range and is defined solely by the calculation convenience. 

This independence of results from /3 additionally manifests itself in the fact that Casimir 
calculations do not actually require any knowledge about the physical volume of the lattice and 
the value of lattice step a. We obtain the dimensionless energy in - units. All geometrical sizes 
of the considered bodies are defined in the lattice step units. Therefore, the physical value of 
the lattice step a fixes all sizes and real physical value of the Casimir energy. 

In noncompact electrodynamics, the rotational symmetry is restored automatically at suf- 
ficiently large lattice distances (this was shown, e.g., in j^). Therefore, the continuous limit 
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is only the transition N ^ oo. In other words, we should find that the correct dependence of 
the Casimir energy on the distance between the bodies is restored at large lattice distances. 

This difference between noncompact QED and non-abelian theories can be understood from 
the following simple arguments. We should decrease the lattice step relatively to a typical scale 
of the theory to achieve the continuous limit. The Compton wavelength of glueball is such a 
scale in non-abelian theories. We vary the relation of the lattice step to this characteristic length 
by changing /3. So we can make the lattice step small in comparison with the glueball scale. 
It means also that the correlation length tends to infinity in the lattice units. The situation is 
quite different in non-compact Abelian theory. There are no glueballs in noncompact Abelian 
theory without fermions in the continuous limit. The correlation length is just a lattice artifact 
of discretization in this theory. The scale in this theory is defined by the sizes of the interacting 
bodies. It means that the continuous limit in this theory is archived in the case when the sizes 
of the bodies (in the lattice units) become large and all the lattice sizes tend to infinity. 
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Figure 2: Calculation of (0|T°°(x)|0)i?'^ in case of the periodic boundary conditions, R is the 
size of the space region in the lattice step units. Lattice sizes: 64 x 64 x 64 x R. The line is 
the analytical answer [Ti] . 

The restoration of the continuous limit in the Casimir problem on the lattice is illustrated in 
Fig. 2. In this calculation the Casimir energy of the electro-magnetic field is studied in the case 
of the space-periodic boundary conditions. This problem provides a good test for our method 
because there is the analytical solution in this case [H]. The lattice observable is the vacuum 
energy density calculated by using the renormalization procedure which we discussed above. 

The energy density of the electro-magnetic field in case of the periodic boundary conditions 
is a constant [H] and is proportional to 1/-R^, where R is the size of the space region where 
we confine the fields by the periodic boundary conditions (this statement is obvious from 
dimensional reasons). Therefore, the quantity {0\T^^[x)\0)R^ should be constant for large size 
of the lattice, when the continuous limit is reproduced. 

In Fig. 2 one can see that our lattice calculation of the energy density tends to the analytical 
result at large lattice distances R. If the lattice distance R is about 11-12 lattice units, the 
lattice results are very close to the analytical ones, and we can conclude that the continuous 
limit is achieved just at these lattice distances. 



7 



We test our method on the well-known problem of the Casimir interaction of two parallel 
conducting planes. In Fig. 3 the lattice calculation of the vacuum energy of the electro- 
magnetic field is presented for this model (we study the planes with quite large size due to the 
continuous limit reasons). We consider two parallel planes which are separated by the distance 
D in 2;-direction. The calculations are performed on the anisotropic lattice with the following 
sizes: 32a at x, y and t directions, Ra^ at z direction, a^/a = 1/3. The analytical result for 
this problem is well known {c = h = 1): 

E 



720D 



where S is the area of the planes, D is the distance between the planes. We fit the calculated 
points by the function E = Pi/B? using the Least Squares method. As one can find, the 
dependence of the Casimir energy on the distance between planes is reproduced correctly by 
our lattice method. In order to compare this lattice result with the analytical answer, one 
should note that all geometrical sizes are proportional to the lattice step a, while the energy is 
proportional to . This means that 

■K^N^ 1 1 

The lattice steps in x, y and t directions are equal to a, but the lattice step is equal to aa in 
deformed z direction. The area of the planes is equal to S* = (A^a)^, where is the lattice size 
in X and y directions; the distance between planes is D = aaR. 
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Figure 3: Calculation of the Casimir interaction energy for two parallel conducting planes. R 
is the distance between planes in the lattice steps: D = aaR. The fitting line is found by the 
Least Squares method: y = j^- All calculation are performed on the anisotropic lattice with 
the deformation coefficient a^jaxyt = 1/3. The lattice sizes: 32a x 32a x 32a x Rqz- 

So we have to compare the dimensionless coefficient at in the expression for the Casimir 
energy fll2l) with the coefficient Pi calculated by fitting the numerical data using the function 
E = Pi/R^. The analytical answer for this coefficient is p^'^°-'-y^ = 379 (the parameters of the 
lattice are shown in Fig. 3). In this test problem one can find that the numerical accuracy of 
our approach is about 5%. 

Let us apply our method to the calculation of the Casimir forces between two conducting 
surfaces with complicated shapes. We consider two rectangular gratings (see Fig. 4). 
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Figure 4: Geometrical parameters of the gratings. 



We perform the calculation for the ideal conducting gratings with the following parameters: 
the width of the "tooth" and the distance between them are Li = L2 = 7a; the height of the 
"tooth" is H = 7a^; the distance between the gratings (between the tops of the "teeth") is 
R = 8a^. All sizes are chosen sufficiently large to restore the continuous limit. The lattice is 
the same that in the previous calculation for two planes: a^/a = 1/3. 

The result of the Casimir energy calculation for the two gratings is shown in Fig. 5. 

We also compare our calculation with the results of approximate methods (PFA and "Op- 
timal" PFA [15]). The PFA approach is often considered as the initial approximation for such 
systems. In the case of the two gratings interaction the PFA method yields the "saw" consisted 
of the straight lines. It is obvious that the deviations from the PFA are rather large, especially 
when the "teeth" of one grating are opposite to the wells of the other one. So a corrections to 
ordinary PFA should be considered. We compare our results with so-called "Optimal" PFA. It 
was proposed by Jaffe and Scardicchio in the paper [15]. In this approach the Casimir energy 
is given by: 

Here D is space region between the boundaries and /i2(r) is the shortest path that passes 
through the point f from one boundary to another. Unfortunately in this case of rectangular 
gratings "Optimal" PFA is not fully adequate too (see Fig. 5). 

We also calculate the vacuum energy density that corresponds to the gratings interaction. 
It is shown in Fig. 6 for different lateral shifts for Monte-Carlo calculations and for "Optimal" 
PFA. 

The calculation of the energy density (see Fig. 6) allows us to clear up the reasons why the 
PFA is not appropriate for this shape of interacting surfaces. In the PFA approach, we take 
into account only interaction of the parts of the surfaces that lie opposite to each other. One 
can see in Fig. 6 that in the situation when the "teeth" are opposite to the "wells" , the vacuum 
energy density is concentrated in the space region between the corners of the "teeth" . And it 
is just the very case when the PFA gives maximal error. Hence we can conclude that the main 
reason for the errors in the PFA is the neglecting of those specific Casimir forces that appear 
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Figure 5: Vacuum energy dependence on the lateral shift of the gratings. Comparison of our 
Monte-Carlo calculations with the results of PFA and "Optimal" PEA methods. Lattice sizes: 
42a X 42a x 42a x 22az. Deformation coefficient: a^ldxyt = 1/3. 



between the "corners" of the gratings. 

"Optimal" PFA takes into account the "corners" interaction. However its weak point in this 
geometry is the edge effect: energy density in the region between parallel parts of the gratings 
strongly depends on the size of this region (see the first column in Fig. 6), but in PFA and 
"Optimal" PFA (see the second column in Fig. 6) energy density in this region is constant. 

Our method can also be considered as Monte-Carlo calculation of some elements of the 
Green's function. So for abelian fields our results will be close to the results obtained by the 
Euclidean Green's function method [3]. 



Conclusions 

We introduced the new method for calculation of the Casimir forces based on the Monte-Carlo 
simulations in noncompact QED on the lattice. The method consists in two parts: the definition 
of the boundary conditions for the field configurations; the definition of the lattice variable that 
gives the vacuum energy density after the averaging over these configurations. 

The method is tested on the problem of the Casimir interaction of two parallel conducting 
plates. Also this approach is used to study the lateral Casimir forces between two rectangular 
gratings. The Casimir energy and the vacuum energy density for the interaction are obtained 
for different lateral shifts of the gratings. 
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Figure 6: Vacuum energy density distribution for different lateral shifts. Results of Monte-Carlo 
calculations (left column) and "Optimal" PFA calculations (right column). 

Acknowledgements 

The authors thank 1. G. Pirozhenko and V. N. Marachevsky for the substantial discussions. 
Computing facilities of the MSU Supercomputer Center have been used for our calculations. 

References 

[1] M. Bordag, G. L. Klimchitskaya, U. Mohideen and V. M. Mostepanenko. Advances in the 
Casimir Effect. Oxford University Press, International Series of Monographs on Physics 
145, 2009. 

[2] M. Bordag and K. Scharnhorst, Phys. Rev. Lett. 81 (1998) 3815. 



11 



[3] M. T. Homer Reid, A. W. Rodriguez, J. White, and S. G. Johnson, Phys. Rev. Lett. 
103 (2009) 040401; M. Levin, A. P. McCauley, A. W. Rodriguez, M. T. Homer Reid, 
and S. G. Johnson, Phys. Rev. Lett. 105 (2010) 090403. 

[4] H. Gies, K. Langfeld, L. Moyaerts, JHEP 018 (2003) 0306; H. Gies and K. Khngmuller, 
Phys. Rev. D 74 (2006) 045002. 

[5] O. Pavlovsky and M. Ulybyshev, Int. J. Mod. Phys. A 25 (2010) 2457 |arXiv:0911.2"635] 
[hep-lat]]. 

[6] O. Pavlovsky, M. Ulybyshev, Particles and Nuclei, Letters, v. 7 (2010) 483. 

[7] A. Lambrecht, V. Marachevsky, Phys. Rev. Lett. 101 (2008) 160403. 

[8] R. Buscher and T. Emig, Phys. Rev. A 69 (2004) 062101. 

[9] D. M. Ceperley, Rev. Mod. Phys. 67 (1995) 279. 
[10] M. Creutz, B. Freedman, Annals of Physics 132 (1981) 427. 
[11] G. Burgers and F. Karsch Nucl. Phys. B 304 (1988) 587. 
[12] T. A. DeGrand and D. Toussaint, Phys. Rev. D 24 (1981) 466. 
[13] Yu. Makeenko, Sov. Phys. Usp. 27 (1984) 401. 
[14] S. G. Mamayev, N. N. Trunov, Theor. Math. Phys. 38 (1979) 228. 
[15] R. L. Jaffe, A. Scardicchio, Phys. Rev. Lett. 92 (2004) 070402. 



12 



